Itinerant ferromagnetism of a repulsive atomic Fermi gas: 
a quantum Monte Carlo study 
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We investigate the phase diagram of a two-component repulsive Fermi gas at T = by means of 
quantum Monte Carlo simulations. For a given value of the positive s-wave scattering length, both 
purely repulsive and purely attractive model potentials are considered in order to analyze the limits 
of the universal regime where the details of interatomic forces can be neglected. The equation of 
state of both balanced and unbalanced systems is calculated as a function of the interaction strength 
and the critical density for the onset of ferromagnetism is determined. The energy per particle of the 
strongly polarized gas is calculated and parametrized in terms of the physical properties of repulsive 
polarons, which are relevant for the stability of the fully magnetized ferromagnetic state. Finally, 
we analyze the phase diagram in the polarization/interaction plane under the assumption that only 
phases with homogeneous magnetization can be produced. 
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Over the past decade there has been substantial 
progress in the experimental realization of quantum de- 
generate atomic Fermi gases. A major part of the activ- 
ity carried out so far was devoted to the investigation 
of the role of attractive interactions, with special em- 
phasis on the onset of pairing and superfluidity in the 
vicinity of a Feshbach resonance as well as in the pres- 
ence of spin imbalance More recently attention was 
drawn to repulsive interactions and the onset of mag- 
netic behavior. This topic is particularly important in 
optical lattices because of its connection with the repul- 
sive Hubbard model, a fundamental paradigm of con- 
densed matter physics with still many unanswered ques- 
tions but also for continuous systems where a major 
recent achievement has been the observation of itiner- 
ant ferromagnetism induced by repulsive forces in a two- 
component Fermi gas This experiment realizes the 
Stoner model, a textbook Hamiltonian that aims to de- 
scribe itinerant ferromagnetism in an electron gas with 
screened Coulomb interaction Q. 

On the theoretical side there have been a number of 
papers addressing the problem of stability of a repulsive 
two-component Fermi gas [j| and of phase separation in 
harmonic trapped configurations within the local density 
approximation [(|. These studies are based on a simple 
mean-field description of interaction effects that is valid 
to linear order in the scattering length. In homogeneous 
systems at T = they predict a second order phase tran- 
sition to a magnetized state if the interaction strength 
is larger than the critical value ij?a > 7r/2, where a is 
the s-wave scattering length and kp — (S^n) 1 / 3 is the 
Fermi wave vector in terms of the total particle density of 
the gas n = n-\ + n\, . An extension of this approach that 
includes next order corrections to the interaction energy 
was developed in Ref. Q and predicts a smaller value of 
the critical density (kpa > 1.054), as well as a discontin- 
uous jump in the magnetization. Low-energy theories of 
itinerant fermions also predict a first-order transition Q . 
A recent non-perturbative quantum Monte Carlo Calcu- 
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FIG. 1: (color online). Phase diagram of the HS gas in the 
interaction/polarization plane. The green region corresponds 
to the homogeneous phase. The other regions correspond to 
phase separated states with partially polarized domains (yel- 
low) and fully ferromagnetic domains (pink). The (blue) sym- 
bols correspond to the minimum of the curve E(P) and the 
solid/dashed line is the phase boundary determined from the 
equilibrium condition for pressure and chemical potentials. 
The blue and red arrows indicate the critical densities where 
X diverges and full ferromagnetism sets in, respectively for 
the HS and SW potential. 



lation, instead, suggests the existence of a textured mag- 
netic phase at the border of the ferromagnetic transition 
and yields the value kpa ~ 0.8 for the critical density 0. 
On the other hand, the existence of a ferromagnetic tran- 
sition has been questioned in Ref. fioj by arguing that 
nonmagnetic states with strong short-ranged repulsive 
correlations could be energetically favorable compared to 
ferromagnetic ones. 

Various important issues concerning the regime of 
strong repulsion are still open. In this Letter we provide 
answer to some of them, in particular: i) we calculate the 
equation of state of the Fermi gas using different poten- 
tials to determine the regime of interaction strength kpa 
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where universality in terms of just the s-wave scattering 
length a is lost and other parameters of the interatomic 
potential become relevant; ii) we study population imbal- 
anced configurations and show that the equation of state 
of the highly polarized gas can be described in terms of a 
Fermi liquid of polarons, similarly to the well-established 
case of attractive interactions [ll|, ; hi) we investigate 
the onset of ferromagnctism and show how the physics of 
polarons is related to the stability of the ferromagnetic 
states; iv) we establish the phase diagram in the polariza- 
tion/interaction plane where we determine the borders of 
three phases fl3j : a uniformly polarized phase and mixed 
phases consisting of partially or fully polarized domains, 
as shown in Fig. [TJ 

To address these issues we perform quantum Monte 
Carlo (QMC) simulations of Fermi systems characterized 
by a positive scattering length a > 0, interacting through 
either purely repulsive or purely attractive forces. The 
Hamiltonian of the Fermi gas is 

ft 2 ( N t N i \ 
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where m denotes the mass of the atoms, 
and label, respectively, spin-up and spin-down 

fcrmions with Nf + N± = N, N being the total number 
of atoms. We model the interspecies interatomic inter- 
actions using three different potentials: i) a hard sphere 
(HS) potential, V{r) = +oo if r < a and zero otherwise, 
ii) a repulsive soft sphere (SS) potential, V(r) = Vq if 
r < Rq and zero otherwise and hi) an attractive square 
well (SW) potential, V(r) = — Vq if r < Rq and zero 
otherwise (Vo > 0). The s-wave scattering length a co- 
incides with the range of the potential in the HS case 
and can readily be determined from the range Rq and 
the strength V in the other two cases. For the SS po- 
tential a is always smaller than the range Rq and we hx 
the height V) such that a = Rq/2. In the case of the 
SW potential instead, the scattering length diverges to 
±oo every time a new bound state enters the well: we hx 
the range such that tiRq = 10 -6 in terms of the particle 
density n and the depth Vo takes values corresponding to 
the positive branch of a with a single bound state. The 
short-range SW potential provides a realistic description 
of interatomic forces in ultracold atoms. 

In the case of purely repulsive interactions we use the 
hxed-node diffusion Monte Carlo (FN-DMC) method. 
This variational method yields an upper bound for the 
ground-state energy of the gas, sampling the lowest- 
energy wave function whose many-body nodal surface 
is the same as that of a trial wave function ipx- FN- 
DMC can give the exact ground-state energy, provided 
one knows the exact nodal surface, and in general the en- 
ergies have been found to be highly accurate even if nodes 
are only approximate (for more details see e.g. [3). Our 
trial wave function is of the Jastrow- Slater form 

v^t(r) = n /(*•«' w#t)jw) . (2) 



where R = (ri, rjv) is the spatial configuration vector 
of the N particles and -D-|-(j,) denotes the Slater determi- 
nant of plane waves in a cubic box of size L with periodic 
boundary conditions, accommodating the iV-f-m particles 
with up (down) spin. The Jastrow correlation term f{r) 
is obtained from the solution of the two-body scattering 
problem with the potential V(r), satisfying the bound- 
ary condition on its derivative f'(r = L/2) = 0. Since 
f(r) > 0, the many-body nodal surface results only from 
the antisymmetric character of ipT and coincides with 
that of a non interacting gas, thus correctly reproducing 
this limit. Another advantage is that the trial function 
([2]) can be used to simulate both unpolarized (iV-j- = ATj.) 
and polarized (JVf > N{) configurations. 

Attractive interactions, as modeled by the SW po- 
tential, are more delicate, because of the presence of 
bound states (molecules) in the true ground state. The 
atomic Fermi gas of interest here is a meta-stable state 
consisting of unbound fcrmionic atoms and no dimcrs 
or other bound molecules. Since we are interested in 
a metastable excited state we cannot use the FN-DMC 
method and resort to the variational Monte Carlo (VMC) 
calculation that provides a stable estimate of the energy 
Evmc = (iPt\H\i(>t}/('4>t\iI)t}- The absence of molecular 
bound states can be readily implemented for two parti- 
cles by choosing the Jastrow correlation term f(r) to be 
the scattering solution of the SW potential corresponding 
to positive energy, which by construction is orthogonal to 
the bound molecule. A many-body wave function is then 
constructed using Eq. ([2]) with this choice of f(r), while 
f'(r = L/2) = takes care of periodic boundary con- 
ditions similarly to the purely repulsive case. For small 
scattering energies, corresponding to large values of L, 
the Jastrow term / changes sign at r = a. The larger 
the size of the simulation box the smaller is the overlap 
between / and the bound-state wave function, and i/jt 
provides an accurate description of the gas-like state in 
the very dilute regime, but in general exhibits a non-zero 
overlap with the state where dimers are formed fill ] . 

It is important to note that the repulsive HS and SS 
potential and the attractive SW potential, even though 
they correspond to the same value of the scattering length 
a and consequently share a similar long-range behavior 
of the correlation functions, exhibit completely different 
short-range correlations as explicitly shown in the inset 
of Fig. [2] where we report results on the antiparallel spin 
pair correlation function. 

Simulations are performed with a maximum number of 
atoms in a single spin component -/V-j-m = 33. For most 
results we checked the influence of finite-size effects by 
repeating the calculations with iVf^ =81. For the re- 
sults of the SW potential, no appreciable change is found 
by reducing the potential range parameter uRq. 

The equation of state of the repulsive dilute gas of 
fcrmions with balanced populations, Nf = iVj,, has pre- 
viously been calculated using perturbation theory to sec- 
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FIG. 2: (color online). Equation of state of the unpolarized 
gas. We show FN-DMC results for the repulsive HS and SS 
potentials and VMC results for the attractive SW potential. 
The perturbation expansion Eq. © is shown with the (black) 
dashed line, while the (green) horizontal line corresponds to 
the energy of the fully ferromagnetic state. Inset: Pair corre- 
lation function g<ti(r). The range of the SS and SW potential 
is respectively Ro = 2a and Ra = 0.06a. The minimum at 
r = a for the SW potential corresponds to the node in the 
Jastrow correlation term f(r). 
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where Ep = Ti 2 kp/2m is the Fermi energy. One should 
notice that higher order terms in the above equation will 
depend not only on the scattering length a, but also on 
other details of the interatomic potential. As can be seen 
in Fig. [2j agreement between the QMC results and Eq. 
(|3"| is found for kpa < 0.4, but significant deviations and 
a gradual loss of universality become evident for larger 
values of the interaction strength. In the figure we also 
show the energy of the fully ferromagnetic (FF) state 
E FF = 3/5E F 2 2 / 3 (Nf + Ni) consisting of two spatially 
separated regions of non-interacting spin-up and spin- 
down fermions. Any uniform mixture of the two spin 
species whose energy exceeds the value E-p-p is clearly 
unstable against phase separation. 

In order to better characterize the critical behavior at 
the onset of ferromagnetic behavior, we calculate the 
equation of state of the gas as a function of the sys- 
tem polarization P = (N t - Ni)/(N f + N±) and then 
show that for P < 1 it can be well described in terms 
of weakly interacting polarons. Results for the HS po- 
tential arc shown in Fig. [3j and similar calculations were 
carried out also for the SW potential (not shown). Fi- 
nite size effects are reduced by subtracting from the QMC 
results, the finite size corrections to the ground state en- 
ergy JSo(JVf, JVj.) — Eq L (P) of non- interacting fermions 
with the same number of particles and the same polar- 
ization P (TL refers here to the thermodynamic limit). 
The validity of this method, that relies on Fermi liquid 



FIG. 3: Energy of the HS function of polarization 

for different values of kpa. The symbols correspond to FN- 
DMC results, while the lines at small and large P correspond, 
respectively, to the fitted quadratic law and polaron energy 
functional of Eq. (j4| . The horizontal dotted line is the thresh- 
old energy Eff of the fully ferromagnetic state. 



theory, is discussed in Ref. [17[ where it is compared with 
twist-averaged boundary conditions. 

From the results at small polarization we extract the 
inverse magnetic susceptibility 1/x using a quadratic fit: 
E(P) = E(P = 0) + NE F (xo/x)P 2 /3, where X o = 
3n/2Ep is the susceptibility of the non-interacting gas. 
The results for both the HS and SW potential are shown 
in the inset of Fig. 2] For kpa > 0.6 we find a linear 
decrease of the inverse magnetic susceptibility with in- 
creasing interaction strength, showing that x diverges at 
the critical densities kpa ~ 0.82 and kpa ~ 0.86, re- 
spectively for the HS and the SW potential. As seen in 
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FIG. 4: (color online). Chemical potential at zero concen- 
tration of the repulsive polaron. We show FN-DMC results 
for the repulsive HS and SS potentials and VMC results for 
the attractive SW potential. The (blue) solid line is a fit to 
the HS results. The second-order perturbation expansion is 
shown with the (black) dashed line, while the (green) hori- 
zontal line corresponds to Ef<\- Inset: Inverse magnetic sus- 
ceptibility 1/x for the HS (blue circles) and SW (red squares) 
potential. The lines are linear fits to the data. 
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Fig. [3] a minimum in E(P) at finite polarization P ap- 
pears when kpa > 0.82. This minimum corresponds to 
a thermodynamically stable mixed state with partially 
polarized domains (see Fig. [TJ. For P < P the system 
follows the coexistence line where, in the thermodynamic 
limit, E(P) is the sum of energies of domains with po- 
larization P and — P whose relative volume changes with 
P. For P > P, on the other hand, a uniform mixture of 
spin components can still survive. 

For large polarizations P the QMC results are well 
fitted by the energy functional of weakly interacting po- 
larons 

E(x) = ^N t E Ft (l + Ax + £^ 5/3 + Fx 2 ^ , (4) 

which results from an expansion in the small concen- 
tration x = N^/Nf of spin-down impurities in a Fermi 
sea of spin-up particles. Here Ep-f = K kp^./2m is the 
Fermi energy of the spin-up particles in terms of the cor- 
responding wave vector kp^ = (Qir 2 ^) 1 ^ . Furthermore, 
the quantities A, m* and F, which are all functions of 
kp^a, denote respectively the chemical potential at zero 
concentration, the effective mass and the interaction pa- 
rameter of the polaronic quasiparticles (see Ref. [ll[ for 
the similar treatment of attractive polarons) . 

The results for A as a function of the interaction 
strength are shown in Fig. [TJ The weak coupling re- 
sult A = Ep^(4:kF^a/3TT + 2(kpfa) 2 /it 2 ), calculated to 
second order in the scattering length [lH, [l9[, is also 
shown. This quantity is important in order to under- 
stand the stability of the fully ferromagnetic phase: in 
fact, if A < Ep^, the fully polarized domain is unstable 
against mixing of spin-up and spin-down particles result- 
ing in a partially polarized state. This problem was in- 
vestigated using a single particle-hole wave function in 
Ref. Jljj . We estimate that the FF state becomes unsta- 
ble at kFfa = 1.2 and kFfO, = 1.4, respectively for the 
HS and SW potential. By imposing pressure and chemi- 



cal potential equilibrium between highly polarized states 
described by the energy functional UJ). we determine the 
critical polarization boundaries, shown in the phase dia- 
gram of Fig. [TJ separating the homogeneous unbalanced 
mixture from the partially ferromagnetic state comprised 
of domains with opposite critical polarization ±P. This 
determination of the phase boundaries is reliable at large 
total polarization P, where the energy functional (UJ) is 
valid, but becomes less accurate at intermediate values of 
P (dashed line in Fig. [TJ and can not be applied at small 
polarization. In the regime of large and intermediate P 
we find good agreement with the critical polarization as 
determined from the minimum of the E{P) curve. The 
critical polarization has a sharp drop to P = close to 
the density where x diverges, not allowing for a clear 
distinction of the order of the phase transition. The pos- 
sible emergence of new phases, such as the spin textured 
phase proposed in Ref. Q, also requires more detailed 
investigations. 

In conclusion we have calculated the equation of state 
of a repulsive gas of fermions as a function of interaction 
strength and spin polarization. We have determined the 
critical density for the onset of ferromagnetic behavior 
and we have investigated the stability of the fully ferro- 
magnetic state. We find, however, that while the qualita- 
tive features of the phase diagrams are well described by 
just the long-range repulsive correlations induced by the 
positive s-wave scattering length, the quantitative deter- 
mination of the phase diagram depends on the details of 
the interaction potential. 

We acknowledge useful discussions with A. Rccati. 
This work was supported by the Swiss National Sci- 
ence Foundation and by a grant from the Army Re- 
search Office with funding from the DARPA OLE pro- 
gram. As part of the European Science Foundation EU- 
ROCORES Program "EuroQUAM-FerMix" it was sup- 
ported by funds from the CNR and the EC Sixth Frame- 
work Programme. 



[1] S. Giorgini, L.P. Pitaevskii and S. Stringari, Rev. Mod. 
Phys. 80, 1215 (2008). 

[2] A. Georges in Ultracold Fermi Gases, Proceedings of the 
Varenna "Enrico Fermi" Summer School, edited by W. 
Ketterle, M. Inguscio and C. Salomon (Societa Italiana 
di Fisica, Bologna, 2007). 

[3] G.-B. Jo et al, Science 325, 1521 (2009). 

[4] E. Stoner, Philos. Mag. 15, 1018 (1933). 

[5] M. Houbiers et al, Phys. Rev. A 56, 4864 (1997). 

[6] M. Amoruso, I. Meccoli, A. Minguzzi and M. Tosi, Eur. 
Phys. J. D 8, 361 (2000); L. Salasnich, B. Pozzi, A. Parola 
and L. Reatto, J. Phys. At. Mol. Opt. Phys. 33, 3943 
(2000); T. Sogo and H. Yabu, Phys. Rev. A 66, 043611 
(2002); J. LeBlanc et al, Phys. Rev. A 80, 013607 (2009). 

[7] R.A. Duine and A.H. MacDonald, Phys. Rev. Lett. 95, 
230403 (2005). 

[8] D. Belitz, T.R. Kirkpatrick and T. Vojta, Phys. Rev. 



Lett. 82, 4707 (1999); D.L. Maslov and A.V. Chubukov, 

Phys. Rev. B 79, 075112 (2009). 
[9] G.J. Conduit, A.G. Green and B.D. Simons, Phys. Rev. 

Lett. 103, 207201 (2009). 
[10] H. Zhai, Phys. Rev. A 80, 051605(R) (2009). 
[11] C. Lobo, A. Recati, S. Giorgini and S. Stringari, Phys. 

Rev. Lett. 97, 200403 (2006). 
[12] Y. Shin et al, Phys. Rev. Lett. 97, 030401 (2006); N. 

Prokof'ev and B. Svistunov, Phys. Rev. B 77, 020408 

(2008); A. Schirotzek, C.-H. Wu, A. Sommer and M.W. 

Zwierlein, Phys. Rev. Le tt. 102, 230402 (2009); S. 

Nascimbene et al, preprint arXiv:0911.0747 
[13] Our analysis has been limited to uniform magnetic 

phases. 

[14] P.J. Reynolds, D.M. Ceperley, B.J. Alder and W.A. 

Lester Jr., J. Chem. Phys. 77, 5593 (1982). 
[15] The overlap with states where three or more-body bound 



5 



states are formed is suppressed if the range Ro of the SW 

potential is small enough. 
[16] K. Huang and C.N. Yang, Phys. Rev. 105, 767 (1957); 

T.D. Lee and C.N. Yang, Phys. Rev. 105, 1119 (1957). 
[17] C. Lin, F.H. Zong and D.M. Ceperley, Phys. Rev. E 64, 



016702 (2001). 

[18] A. Recati, private communication; C. Mora and F. 

Chevy, preprint [arXiv:1003.02i3l 

[19] X. Cui and H. Zhai, preprint larXiv:1001. 51391 



